########################################################################
################ R CODE FOR ABBOTT, REICHERT and SASMAZ ################
########################################################################

library("foreign")
library("sandwich")
library("stargazer")
library("lmtest")
library("mvtnorm")
library("boot")
library("Zelig")
library("ZeligChoice")


setwd("C:/Users/Jared/SkyDrive/Harvard/Gov 2001/paper")



load("data.const.wnewvar.Rda")


#######################################################################
# Descriptive variables table in Data & Empirical Strategy section    #
#######################################################################

## Table 1 
stargazer(data.const.wnewvar[c("share1831", "whig_change", "seat_change", 
                      "constriots10",
                      "share1826", "share1826_2", "reform_1830",
                      "county", "university", "narrow", "rottenindex",
                      "h_index1831", "primary", "secondary", "tetitary",
                      "population1831", "density1831", "economictrend1", "economictrend2",
                      "contested1", "prop_value_percap")], 
          covariate.labels=c("Whig share 1831 (%)", "Vote change for Whigs", "Seat change for Whigs",
                             "Riots within 10 km",
                             "Whig share 1826", 
                             "Whig share 1826 squared", "Reform support 1830",
                             "County constituency", "University constituency",
                             "Narrow franchise", "Patronage index",
                             "Emp. fract. index", "Agriculture (emp.share)",
                             "Trade (emp.share)", "Professionals (emp. share)",
                             "Population (in 1000s)", "Population density", 
                             "Thriving economy", "Declining economy",
                             "Contested elections", "Income level per capita"), 
          digits=2)


#######################################################################################
##     Tables in the Results Section                                                 ##
#######################################################################################

############### Table 2 #################

model_1.1 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar)
summary(model_1.1)

model_1.2 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap + contested1 + constriots10*contested1, 
                data=data.const.wnewvar)
summary(model_1.2)

model_1.3 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar[ which(data.const.wnewvar$contested1 == 0),])
summary(model_1.3)

model_1.4 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar[ which(data.const.wnewvar$contested1 == 1),])
summary(model_1.4)

stargazer(model_1.1, model_1.2, model_1.3, model_1.4, type = "text",  
          dep.var.labels="Whig Share 1831 (%)",
          covariate.labels=c("Riots within 10 km",
                             "Whig share 1826", 
                             "Whig share 1826 squared", "Reform support 1830",
                             "County constituency", "University constituency",
                             "Narrow franchise", "Patronage index",
                             "Emp. fract. index", "Agriculture (emp.share)",
                             "Trade (emp.share)", "Professionals (emp. share)",
                             "Population (in 1000s)", "Population density", 
                             "Thriving economy", "Declining economy", "Income level per capita", "Contested districts", 
                             "Riots within 10km * Contested dist."), 
          digits=2)



############### Table 3 ##################

model_2.1 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar)
summary(model_2.1)

model_2.2 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap + contested1 + constriots10*contested1, 
                data=data.const.wnewvar)
summary(model_2.2)

model_2.3 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar[ which(data.const.wnewvar$contested1 == 0),])
summary(model_2.3)

model_2.4 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                  county + university + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar[ which(data.const.wnewvar$contested1 == 1),])
summary(model_2.4)

stargazer(model_2.1, model_2.2, model_2.3, model_2.4, type = "text",
          dep.var.labels="Seat change for Whigs",
          covariate.labels=c("Riots within 10 km",
                             "Whig share 1826", 
                             "Whig share 1826 squared", 
                             "County constituency", "University constituency",
                             "Narrow franchise", "Patronage index",
                             "Emp. fract. index", "Agriculture (emp.share)",
                             "Trade (emp.share)", "Professionals (emp. share)",
                             "Population (in 1000s)", "Population density", 
                             "Thriving economy", "Declining economy", "Income level per capita", "Contested districts", 
                             "Riots within 10km * Contested dist."),
          digits=3)



############### Table 4 ##################

model_3.1 <- lm(whig_change ~ constriots10 + share1826 + share1826_2 +
                  county + narrow + rottenindex + h_index1831 + 
                  primary + secondary + tetitary + population1831 + density1831 +
                  economictrend1 + economictrend2 + prop_value_percap, 
                data=data.const.wnewvar[ which(data.const.wnewvar$contested1 == 1),])
summary(model_3.1)

stargazer(model_3.1, type = "text",
          dep.var.labels="Vote change for Whigs",
          covariate.labels=c("Riots within 10 km",
                             "Whig share 1826", 
                             "Whig share 1826 squared", "County constituency", 
                             "Narrow franchise", "Patronage index",
                             "Emp. fract. index", "Agriculture (emp.share)",
                             "Trade (emp.share)", "Professionals (emp. share)",
                             "Population (in 1000s)", "Population density", 
                             "Thriving economy", "Declining economy",
                             "Income level per capita"), 
          digits=3)



###########################################################################
###################### SPATIAL VARIABLES SECTION ##########################
###########################################################################

SpatialVar2 <- read.csv("SpatialVar2.csv")
data.const.wnewvar4 <- merge(data.const.wnewvar, SpatialVar2, by = "consid")

## Descriptive variables table in the Spatiality section - Table 5
stargazer(data.const.wnewvar4[c("near_comp_10", "near_comp_20", "near_comp_30", "near_comp_40", "Distance")], 
          type = "text",
          covariate.labels=c("Contested districts within 10 km",
                             "Contested districts within 20 km",
                             "Contested districts within 30 km",
                             "Contested districts within 40 km",
                             "Distance to the nearest contested district"), 
          digits=2)


#### Spatial Variable Regressions - Table 6
data.const.wnewvar4 <- merge(data.const.wnewvar, SpatialVar2, by = "consid")
f.2a5s <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
               county + university + narrow + rottenindex + h_index1831 + 
               primary + secondary + tetitary + population1831 + density1831 +
               economictrend1 + economictrend2, data=data.const.wnewvar4)


f.2a5s1 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_10, data=data.const.wnewvar4)

f.2a5s2 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_20, data=data.const.wnewvar4)

f.2a5s3 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_30, data=data.const.wnewvar4)

f.2a5s4 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_40, data=data.const.wnewvar4)

f.2a5s7 <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + Distance, data=data.const.wnewvar4)


stargazer(f.2a5s, f.2a5s1, f.2a5s2, f.2a5s3, f.2a5s4, f.2a5s7, 
          dep.var.labels="Whig Share 1831 (%)", font.size="tiny",
          covariate.labels=c("Riots within 10 km", "Whig share 1826", 
                             "Whig share 1826 squared", "Reform support 1830",
                             "County constituency", "University constituency", 
                             "Narrow franchise", "Patronage index",
                             "Emp. fract. index", "Agriculture (emp.share)", 
                             "Trade (emp. share)", "Professionals (emp. share)",
                             "Population", "Population density", "Thriving economy", 
                             "Declining economy", "Contested dist. within 10 km", "Contested dist. within 20 km",
                             "Contested dist. within 30 km", "Contested dist. within 40 km", "Distance to the nearest contested dist.")
)

## Quantity of interest with Zelig
fz.2a5s4 <- zelig(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                    county + university + narrow + rottenindex + h_index1831 + 
                    primary + secondary + tetitary + population1831 + density1831 +
                    economictrend1 + economictrend2 + near_comp_40, model="ls", data=data.const.wnewvar4)

xlow <- setx(fz.2a5s4, near_comp_40=0)
xhigh <- setx(fz.2a5s4, near_comp_40=2)
sim.2a5s4 <- sim(fz.2a5s4, x=xlow, x1=xhigh)
summary(sim.2a5s4)
plot(sim.2a5s4)

#################################################################
####### APPENDIX ################################################
#################################################################

#### Ordered probit - Table A1 

fit.2a6op <- zelig(
  factor(share1831) ~ constriots10 + share1826 + reform_1830 + county + university + rottenindex + 
    economictrend2, 
  data=data.const.wnewvar[which(data.const.wnewvar$share1831 != 25 & data.const.wnewvar$share1831 != 75),], 
  model="oprobit")

fit.2a6opwci <- zelig(
  factor(share1831) ~ constriots10 + share1826 + reform_1830 + county + university + rottenindex + 
    economictrend2 + prop_value_percap + contested1 + constriots10*contested1, 
  data=data.const.wnewvar[which(data.const.wnewvar$share1831 != 25 & data.const.wnewvar$share1831 != 75),], 
  model="oprobit")

fit.2a6opwcii <- zelig(
  factor(share1831) ~ constriots10 + share1826 + reform_1830 + county + university + rottenindex + 
    economictrend2 + prop_value_percap, 
  data=data.const.wnewvar[which(data.const.wnewvar$share1831 != 25 & data.const.wnewvar$share1831 != 75 & data.const.wnewvar$contested1 == 0),], 
  model="oprobit")

fit.2a6opwciii <- zelig(
  factor(share1831) ~ constriots10 + share1826 + reform_1830 + county + university + rottenindex + 
    economictrend2 + prop_value_percap, 
  data=data.const.wnewvar[which(data.const.wnewvar$share1831 != 25 & data.const.wnewvar$share1831 != 75 & data.const.wnewvar$contested1 == 1),], 
  model="oprobit")

stargazer(fit.2a6op, fit.2a6opwci, fit.2a6opwcii, fit.2a6opwciii,
          dep.var.labels="Whig Share 1831 (%)",
          covariate.labels=c("Riots within 10 km",
                             "Whig share 1826", 
                             "Whig share 1826 squared", "Reform support 1830",
                             "County constituency", "University constituency",
                             "Patronage index",
                             "Declining economy", "Income level per capita", "Contested districts", 
                             "Riots within 10km * Contested dist."))

#### Data validation - Table A2 and corresponding figures (Figure A1, A2 and A3)


data.const.wnewvar.trim <- data.const.wnewvar[is.na(data.const.wnewvar$prop_value_percap) == F,]
data.const.wnewvar.trim.c <- data.const.wnewvar.trim[data.const.wnewvar.trim$contested1 == 1,]
data.const.wnewvar.trim.uc <- data.const.wnewvar.trim[data.const.wnewvar.trim$contested1 == 0,]

###################### models from table 2 ##################################

EVerr.m1.1 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim[sample(nrow(data.const.wnewvar.trim)),]
  train.data <- data_shuffled[1:166,]
  test.data <- data_shuffled[166:208,]
  
  model <- lm(share1831 ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m1.1[j] <- mean(errors)            
}
SEmin.m1.1 <- mean(EVerr.m1.1) - ((0.95 * (max(EVerr.m1.1) - min(EVerr.m1.1)))/2)
SEmax.m1.1 <- mean(EVerr.m1.1) + ((0.95 * (max(EVerr.m1.1) - min(EVerr.m1.1)))/2)

EVerr.m1.2 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim[sample(nrow(data.const.wnewvar.trim)),]
  train.data <- data_shuffled[1:166,]
  test.data <- data_shuffled[166:208,]
  
  model <- lm(share1831 ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap + constriots10*contested1, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m1.2[j] <- mean(errors)            
}
SEmin.m1.2 <- mean(EVerr.m1.2) - ((0.95 * (max(EVerr.m1.2) - min(EVerr.m1.2)))/2)
SEmax.m1.2 <- mean(EVerr.m1.2) + ((0.95 * (max(EVerr.m1.2) - min(EVerr.m1.2)))/2)

EVerr.m1.3 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim.uc[sample(nrow(data.const.wnewvar.trim.uc)),]
  train.data <- data_shuffled[1:138,]
  test.data <- data_shuffled[139:173,]
  
  model <- lm(share1831 ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m1.3[j] <- mean(errors)            
}
SEmin.m1.3 <- mean(EVerr.m1.3) - ((0.95 * (max(EVerr.m1.3) - min(EVerr.m1.3)))/2)
SEmax.m1.3 <- mean(EVerr.m1.3) + ((0.95 * (max(EVerr.m1.3) - min(EVerr.m1.3)))/2)

EVerr.m1.4 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim.c[sample(nrow(data.const.wnewvar.trim.c)),]
  train.data <- data_shuffled[1:28,]
  test.data <- data_shuffled[29:35,]
  
  model <- lm(share1831 ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m1.4[j] <- mean(errors)            
}
SEmin.m1.4 <- mean(EVerr.m1.4) - ((0.95 * (max(EVerr.m1.4) - min(EVerr.m1.4)))/2)
SEmax.m1.4 <- mean(EVerr.m1.4) + ((0.95 * (max(EVerr.m1.4) - min(EVerr.m1.4)))/2)

###################### models from table 3 ##################################

EVerr.m2.1 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim[sample(nrow(data.const.wnewvar.trim)),]
  train.data <- data_shuffled[1:166,]
  test.data <- data_shuffled[166:208,]
  
  model <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m2.1[j] <- mean(errors)            
}
SEmin.m2.1 <- mean(EVerr.m2.1) - ((0.95 * (max(EVerr.m2.1) - min(EVerr.m2.1)))/2)
SEmax.m2.1 <- mean(EVerr.m2.1) + ((0.95 * (max(EVerr.m2.1) - min(EVerr.m2.1)))/2)

EVerr.m2.2 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim[sample(nrow(data.const.wnewvar.trim)),]
  train.data <- data_shuffled[1:166,]
  test.data <- data_shuffled[166:208,]
  
  model <- lm(seat_change ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap + constriots10*contested1, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m2.2[j] <- mean(errors)            
}
SEmin.m2.2 <- mean(EVerr.m2.2) - ((0.95 * (max(EVerr.m2.2) - min(EVerr.m2.2)))/2)
SEmax.m2.2 <- mean(EVerr.m2.2) + ((0.95 * (max(EVerr.m2.2) - min(EVerr.m2.2)))/2)

EVerr.m2.3 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim.uc[sample(nrow(data.const.wnewvar.trim.uc)),]
  train.data <- data_shuffled[1:138,]
  test.data <- data_shuffled[139:173,]
  
  model <- lm(seat_change ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m2.3[j] <- mean(errors)            
}
SEmin.m2.3 <- mean(EVerr.m2.3) - ((0.95 * (max(EVerr.m2.3) - min(EVerr.m2.3)))/2)
SEmax.m2.3 <- mean(EVerr.m2.3) + ((0.95 * (max(EVerr.m2.3) - min(EVerr.m2.3)))/2)

EVerr.m2.4 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim.c[sample(nrow(data.const.wnewvar.trim.c)),]
  train.data <- data_shuffled[1:28,]
  test.data <- data_shuffled[29:35,]
  
  model <- lm(seat_change ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m2.4[j] <- mean(errors)            
}
SEmin.m2.4 <- mean(EVerr.m2.4) - ((0.95 * (max(EVerr.m2.4) - min(EVerr.m2.4)))/2)
SEmax.m2.4 <- mean(EVerr.m2.4) + ((0.95 * (max(EVerr.m2.4) - min(EVerr.m2.4)))/2)


###################### models from table 4 ##################################

EVerr.m3.1 <- NA
for(j in 1:50){
  data_shuffled <- data.const.wnewvar.trim.c[sample(nrow(data.const.wnewvar.trim.c)),]
  train.data <- data_shuffled[1:28,]
  test.data <- data_shuffled[29:35,]
  
  model <- lm(whig_change ~ share1826 + share1826_2 + reform_1830 +
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + prop_value_percap, data = train.data)
  
  pr.values <- rep(NA, nrow(test.data))
  for (i in 1:nrow(test.data)){
    pr.values[i]<- predict(model, newdata = test.data[i,])  }
  
  errors <- pr.values - test.data$share1831
  EVerr.m3.1[j] <- mean(errors)            
}
SEmin.m3.1 <- mean(EVerr.m3.1) - ((0.95 * (max(EVerr.m3.1) - min(EVerr.m3.1)))/2)
SEmax.m3.1 <- mean(EVerr.m3.1) + ((0.95 * (max(EVerr.m3.1) - min(EVerr.m3.1)))/2)




################## figure for models in table 2 ###############
plot(mean(EVerr.m1.1), 1,
     pch = 20,
     col = "dark blue",
     ylim = c(0, 10),
     xlim = c(-100, 25),
     xlab = "mean error (predicted values - actual values)",
     ylab = "model",
     main = "Cross-Validation Results by Model")
lines(x = c(SEmin.m1.1, SEmax.m1.1), y = c(1, 1),
      col = "dark blue")

points(mean(EVerr.m1.2), 2,
       pch = 20,
       col = "dark blue")
lines(x = c(SEmin.m1.2, SEmax.m1.2), y = c(2, 2),
      col = "dark blue")

points(mean(EVerr.m1.3), 3,
       pch = 20,
       col = "dark blue")
lines(x = c(SEmin.m1.3, SEmax.m1.3), y = c(3, 3),
      col = "dark blue")

points(mean(EVerr.m1.4), 4,
       pch = 20,
       col = "dark blue")
lines(x = c(SEmin.m1.4, SEmax.m1.4), y = c(4, 4),
      col = "dark blue")

################## figure for models in table 3 ###############


plot(mean(EVerr.m2.1), 5,
     pch = 20,
     col = "dark red")
lines(x = c(SEmin.m2.1, SEmax.m2.1), y = c(5, 5),
      col = "dark red")

points(mean(EVerr.m2.2), 6,
       pch = 20,
       col = "dark red")
lines(x = c(SEmin.m2.2, SEmax.m2.2), y = c(6, 6),
      col = "dark red")

points(mean(EVerr.m2.3), 7,
       pch = 20,
       col = "dark red")
lines(x = c(SEmin.m2.3, SEmax.m2.3), y = c(7, 7),
      col = "dark red")

points(mean(EVerr.m2.4), 8,
       pch = 20,
       col = "dark red")
lines(x = c(SEmin.m2.4, SEmax.m2.4), y = c(8, 8),
      col = "dark red")

################## figure for models in table 4 ###############

plot(mean(EVerr.m3.1), 9,
     pch = 20,
     col = "goldenrod3")
lines(x = c(SEmin.m3.1, SEmax.m3.1), y = c(9, 9),
      col = "goldenrod3")




#### Spatiality - doing with the seat change variable - Table A3

f.2a5t <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
               county + university + narrow + rottenindex + h_index1831 + 
               primary + secondary + tetitary + population1831 + density1831 +
               economictrend1 + economictrend2, data=data.const.wnewvar4)


f.2a5t1 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_10, data=data.const.wnewvar4)

f.2a5t2 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_20, data=data.const.wnewvar4)

f.2a5t3 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_30, data=data.const.wnewvar4)

f.2a5t4 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + near_comp_40, data=data.const.wnewvar4)

f.2a5t7 <- lm(seat_change ~ constriots10 + share1826 + share1826_2 + 
                county + university + narrow + rottenindex + h_index1831 + 
                primary + secondary + tetitary + population1831 + density1831 +
                economictrend1 + economictrend2 + Distance, data=data.const.wnewvar4)

stargazer(f.2a5t, f.2a5t1, f.2a5t2, f.2a5t3, f.2a5t4, f.2a5t7, 
          dep.var.labels="Seat Change for Whigs", font.size="tiny",
          covariate.labels=c("Riots within 10 km", "Whig share 1826", 
                             "Whig share 1826 squared",
                             "County constituency", "University constituency", 
                             "Narrow franchise", "Patronage index",
                             "Emp. fract. index", "Agriculture (emp.share)", 
                             "Trade (emp. share)", "Professionals (emp. share)",
                             "Population", "Population density", "Thriving economy", 
                             "Declining economy", "Contested dist. within 10 km", "Contested dist. within 20 km",
                             "Contested dist. within 30 km", "Contested dist. within 40 km", "Distance to the nearest contested dist.")
)
